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ABSTRACT 

In preparation for an experimental study of magnetorotational instability 
(MRI) in liquid metal, we present non-ideal two-dimensional magnetohydrody- 
namic simulations of the nonlinear evolution of MRI in the experimental geome- 
try. The simulations adopt initially uniform vertical magnetic fields, conducting 
radial boundaries, and periodic vertical boundary conditions. No-slip conditions 
are imposed at the cylinders. Our linear growth rates compare well with existing 
local and global linear analyses. The MRI saturates nonlinearly with horizon- 
tal magnetic fields comparable to the initial axial field. The rate of angular 
momentum transport increases modestly but significantly over the initial state. 
For modest fluid and magnetic Reynolds numbers Re, Rm ~ 10^ — 10^, the fi- 
nal state is laminar reduced mean shear except near the radial boundaries, and 
with poloidal circulation scaling as the square root of resistivity, in partial agree- 
ment with the analysis of Knobloch and Julien. A sequence of simulations at 
Rm = 20 and 10^ < Re < lO'^'^ enables extrapolation to the experimental regime 
{Rm ~ 20, Re ~ 10^), albeit with unrealistic boundary conditions. MRI should 
increase the experimentally measured torque substantially over its initial purely 
hydrodynamic value. 

Subject headings: accretion, accretion disk — instability — (magnetohydrodynamics:) 
MHD — methods: numerical 
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1. Introduction 

Rapid angular momentum transport in accretion disks has been a longstanding astro- 
physical puzzle. The molecular viscosity of astrophysical gases and plasmas is completely 
inadequate to explain observationally inferred accretion rates, so that a turbulent viscosity 
is required. Recent theoretical work (Pringle 1981; Balbus & Hawley 1991; Hawley & Balbus 
1991; Balbus & Hawley 1998) indicates that purely hydrodynamic instabilities are absent 
or ineffective, but that magnetorotational instabilities (MRI) are robust and support vigor- 
ous turbulence in electrically-conducting disks. Although originally discovered by Vclikhov 
(1959) and Chandrasekhar (1960), MRI did not come to the attention of the astrophysical 
community until rediscovered by Balbus & Hawley (1991) and verified numerically (Hawley 
et al. 1995; Brandenburg et al. 1995; Matsumoto & Tajima 1995). It is now believed that 
MRI drives accretion in disks ranging from quasars and X-ray binaries to cataclysmic vari- 
ables and perhaps even protoplanetary disks (Balbus & Hawley 1998). Some astrophysicists, 
however, argue from laboratory evidence that purely hydrodynamic turbulence may account 
for observed accretion rates, especially in cool, poorly conducting disks where MRI may not 
operate (DubruUe 1993; Richard & Zahn 1999; Duschl et al. 2000; Hure et al. 2001). 

Although its existence and importance are now accepted by most astrophysicists, MRI 
has yet to be clearly demonstrated in the laboratory, notwithstanding the claims of Sisan 
et al. (2004), whose experiment proceeded from a background state that was not in MHD 
equilibrium. Recently (Ji et al. 2001; Goodman & Ji 2002), wc have therefore proposed an 
experimental study of MRI using a magnetized Couette flow: that is, a conducting liquid 
(gallium) bounded by concentric differentially rotating cylinders and subject to an axial 
magnetic field. The radii of the cylinders are ri < r2, as shown in Fig. 1; their angular 
velocities, fli & fl2, have the same sign in all cases of interest to us. If the cylinders 
were infinitely long — very easy to assume theoretically, but rather more difficult to build 
experimentally — the steady-state solution would be ideal Taylor-Couette fiow: 

Q(r)=a + A (1) 

where a = (r22'"2 / i'^2~'^i) ^^'^ ^ ~ ''"i^'il^i ~ ^2) / (^"2 ~ '''i)- unmagnetized and 

inviscid limit, such a flow is linearly axisymmetric stable if and only if the specific angular 
momentum increases outwards: that is, {flirfy < (02?'2)^) or equivalently, a6 > 0. A vertical 
magnetic field may destabilize the fiow, however, provided that the angular velocity decreases 
outward, Q| < Q^; in ideal MHD, instability occurs at arbitrarily weak field strengths (Balbus 
& Hawley 1991). The challenge for experiment, however, is that hquid-metal fiows are very 
far from ideal on laboratory scales. While the fiuid Reynolds number Re = Jliri(r2 —ri)/v 
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can be large, the corresponding magnetic Reynolds number 

Rm ^ ^hn(ll^ (2) 
V 

is modest or small, because the magnetic Prandtl number Pm = u/r] 10^® in liquid metals. 
Standard MRI modes will not grow unless both the rotation period and the Alfven crossing 
time are shorter than the timescale for magnetic diffusion. This requires both Rm > 1 and 
S >!, where 

g^VAr^ (3) 
rj 

is the Lundquist number, and Va = B / ^Anp is the Alfven speed. Therefore, Re > 10^ and 
fields of several kilogauss must be achieved in typical experimental geometries. 

Recently, it has been discovered that MRI modes may grow at much reduced Rm and S 
in the presence of a helical background field, a current-free combination of axial and toroidal 
field (Hollerbach & Riidiger 2005; Riidiger et. al. 2005). We have investigated these helical 
MRI modes. While we confirm the quantitative results given by the authors just cited for the 
onset of instability, we have uncovered other properties of the new modes that cast doubt 
upon both their experimental realizability and their relevance to astrophysical disks. To 
limit the length of the present paper, we present results for purely axial background fields 
only. A paper on helical MRI is in preparation. 

One may question the relevance of experimental to astrophysical MRI, especially its 
nonlinear phases. In accretion disks, differential rotation arises from radial force balance 
between the gravitational attraction of the accreting body and centrifugal force. Thermal and 
magnetic energies are small compared to orbital energies, at least if the disk is vertically thin 
compared to its radius. Consequently, nonlinear saturation of MRI cannot occur by large- 
scale changes in rotation profile. In experiments, however, differential rotation is imposed 
by viscous or other weak forces, and the incompressiblity of the fluid and its conflncmcnt 
by a container allow radial force balance for arbitrary VL{r). Thus, saturation may occur by 
reduction of differential rotation, which is the source of free energy for the instability. In this 
respect, MRI experiments and the simulations of this paper may have closer astrophysical 
counterparts among differentially rotating stars, where rotation is subsonic and boundaries 
are nearly stress- free (Balbus & Hawley 1994; Menou, Balbus, & Spruit 2005). Both in the 
laboratory and in astrophysics, however, nonlinear MRI is expected to enhance the radial 
transport of angular momentum. Quantifying the enhanced transport in a Couette flow is a 
primary goal of the Princeton MRI experiment and of the present paper. 

Another stated goal of the Princeton experiment is to vahdate astrophysical MHD codes 
in a laboratory setting. Probably the most widely used astrophysical MHD code is ZEUS 
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(Stone Sz Norman 1992a,b), which exists in several variants. The simulations of this paper 
use ZEUS-2D. Like most other astrophysical MHD codes, ZEUS-2D was designed for com- 
pressible, ideal-MHD flow with simple boundary conditions: outflow, inflow, reflecting — but 
not no-slip. ZEUS would not be the natural choice of a computational fluid-dynamicist 
interested in Couette flow for its own sake. Nevertheless, after modifying ZEUS-2D to incor- 
porate resistivity, viscosity, and no-slip boundary conditions, we find it to be a robust and 
flexible tool for the subsonic flows of interest to us. It reproduces the growth rates predicted 
for incompressible flow (§3), and agrees with hydrodynamic laboratory data (Burin et al. 
2005); MHD data are not yet available. Of course, all real flows are actually compressible; in 
an ideal gas of fixed total volume, density changes generally scale ~ when Mach number 
^ow/^ound < 1- Incompressibility is an idealization in the limit M — 0. We have used an 
isothermal equation of state in ZEUS with a sound speed chosen so that the maximum of 
M < 1/4 and obtain quantitative agreement with incompressible codes at the few-percent 
level (§2). 

Most of the parameters of the simulations in §§3-4 are chosen to match those of the 
experiment. We adopt the same cylinder radii (Fig. 1). The experimental rotation rates of 
both cylinders (and of the endcaps) are separately adjustable, as is the axial magnetic field. 
For these simulations, we adopt fixed values within the achievable range: fli = 4000 rpm & 
0,2 — 533 rpm, B^o = 5000 G. We set the density of the fluid to that of gallium, p = 6 g cm~^. 

Our simulations depart from experimental reality in two important respects: Reynolds 
number and vertical boundary conditions. Computations at Re > 10^ are out of reach of 
any present-day code and computer, at least in three dimensions; Re ~ 10^ might just be 
achievable in axisymmetry, but higher- it!e flows are more likely to be three-dimensional, so 
that an axisymmetric simulation at such a large Re is of doubtful relevance. (The same 
objection might be leveled at all of our simulations for Re ^ 10^. Those simulations are 
nevertheless useful for establishing scaling relations, even if the applicability of the relations 
to real three-dimensional flows is open to question.) We use an artiflcially large kinematic 
viscosity so that Re = 10^ — 10'*''', whereas for the true kinematic viscosity of gallium 
(z/ ^ 3 X 10^^ cni^s^^). Re ^ 10^ at the dimensions and rotation rates cited above. In 
defense of this approximation, we point to the fact that extrapolations of Ekman-circulation 
rates and rotation proflles simulated at Re < 10^ agree well with measurements taken at 
Re — 10^ both in a prototype experiment (Kageyama et al. 2004), and in the present 
aparatus (Burin et al. 2005). We are able to reproduce the experimental values of the 
dimensionless parameters based on resistivity: Rm ~ 20, 5" ~ 4; we also report simulations 
at Rm ~ 10^ - 10^. (The actual diffusivity of gaUium is 77 ~ 2 x lO^cm^s"-^). 

Except for hydrodynamic test simulations carried out to compare with incompressible 
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results and laboratory data (§2), we adopt vertically periodic boundary conditions for all 
fluid variables, with a periodicity length — 2h, where h = 27.9 cm is the actual height 
of the experimental flow. Such boundary conditions are physically unrealistic, but almost 
all published linear analyses of MRl in Couette flows have adopted them because they per- 
mit a complete separation of variables (Ji et al. 2001; Goodman & Ji 2002; Noguchi et al. 
2002; Riidiger & Shalybkov 2002; Riidiger et al. 2003); an exception is Riidiger & Zhang 
(2001). Thus by adopting periodic vertical boundaries, we are able to test our code against 
well-established linear growth rates and to explore — apparently for the first time in Cou- 
ette geometry — the transition from linear growth to nonlinear saturation. The imposition of 
no-slip conditions at finite endcaps introduces important complications to the basic state, in- 
cluding Ekman circulation and Stewartson layers, which we are currently studying, especially 
as regards their modification by the axial magnetic field. But the experimental apparatus 
has been designed to minimize these complications {e.g. by the use of independently con- 
trolled split endcaps) in order to approximate the idealized Couette fiows presented here, 
whose nonlinear development already presents features of interest. This paper is the first 
in a series; later papers will address the effects of finite endcaps on magnetized flow, helical 
MRI instabihties, etc. 



ZEUS-2D offers the option of cartesian {x,y), spherical {R,9), or cylindrical {z,r) co- 
ordinates. We use {z,r). Although all quantities are assumed independent of the azimuth 
ifi, the azimuthal components of velocity (v^) and magnetic fleld (B^) are represented. We 
have implemented vertically periodic boundary conditions (period^ 2h) for all variables, and 
conducting radial boundary conditions for the magnetic fleld. Impenetrable, no-slip radial 
boundaries are imposed on the velocities. Viscosity and resistivity have been added to the 
code. In order to conserve angular momentum precisely, we cast the azimuthal component 
of the Navier-Stokes equation in conservative form: 



in which L — rV^, and Fr and are the viscous angular- momentum fluxes per unit mass. 



2. 



Modifications to ZEUS-2D and Code Tests 








In the spirit of ZEUS, the viscous part of eq. (4) is implemented as part of the "source" 
substep. In accord with the Constrained Transport algorithm (Evans & Hawley 1988), which 
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preserves V • B = 0, resistivity is implemented by an ohmic term added to the electromotive 
force, which becomes 

£ = V X B - 77V X B . (6) 



2.1. Code Tests (1) - Wendl's Low-Re Solution 

At i?e <S 1 and Rm — 0, poloidal flow is negligible and the toroidal flow is steady. V^p 
satisfles ^ 

^(v' - :;:^)v^ = 0. (7) 

Wendl (1999) has given the analytic solution of this equation for no-slip vertical boundaries 
co-rotating with the outer cylinder. This serves as one benchmark for the viscous part of our 
code; note that the vertical boundary conditions differ from those used in the simulations of 
§3-4. 

Figure 2 compares results from ZEUS-2D with the analytical result. The maximum 
relative error is less than 3%. We have also calculated the viscous torque across the mean 
cylinder (r = (ri + r2)/2). Wendl's solution predicts —1.5004 x lO^gcm^s"^, and our 
simulations yield —1.5028 x lO^gcm^s"^. 



2.2. Code Tests (2) - Magnetic Diffusion 

If the fluid is constrained to be at rest, then the toroidal induction equation becomes 

dt ^ \ dr'^ r dr dz'^ J 

An exact solution compatible with our boundary conditions is: 

B = e^S° + e^-^ cos{kz) exp{-rikH) (9) 

where k is the wave number, and and B^ are constants. 

A comparison of the theoretical and simulated results shows that the error scales 
quadratically with cell size, as expected for our second-order difference scheme (Table 1). 
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2.3. Compcirison with an Incompressible Code 

ZEUS-2D is a compressible code. However our experimental fluid, gallium, is nearly 
incompressible at flow speeds of interest, which are much less than its sound speed, 2.7kms~^. 
As mentioned in §1, we can approximate incompressible flow by using a subsonic Mach 
number, M < 1. However, since ZEUS is explicit, M <^ 1 requires a very small time step 
to satisfy the CFL stability criterion. As a compromise, we have used M — l/A (based on 
the inner cylinder) throughout all the simulations presented in this paper. We assume an 
isothermal equation of state to avoid increases in M by viscous and resistive heating; the 
nonlinear compressibility and thermodynamic properties of the actual liquid are in any case 
very different from those of ideal gases, for which ZEUS was written. Figure 3 compares 
results obtained from ZEUS-2D with simulations performed by Kageyama et al. (2004) using 
their incompressible Navier-Stokes code. 



3. Linear MRI Simulations 

In the linear regime, MRI has been extensively studied both locally and globally (Ji 
et al. 2001; Goodman & Ji 2002; Riidiger & Zhang 2001; Noguchi et al. 2002; Riidiger & 
Shalybkov 2002; Riidiger et al. 2003). We have used these linear results to benchmark our 
code. 

In the linear analyses cited above, the system is assumed to be vertically periodic with 
periodicity length 2h, twice the height height of the cylinders. In cylindrical coordinates, the 
equilibrium states are Bq = B^Bz and Vq = rQe^p. WKB methods describe the stability of 
this system very well even on the largest scales (Ji et al. 2001; Goodman & Ji 2002). Linear 
modes are proportional to exp(7t — ikzz)f{krr), where 7 is the growth rate, and f{x) is an 
approximately sinusoidal radial function, at least outside boundary layers, whose zeros are 
spaced by Ax ~ tt. The wavenumbers — mr/h and kr ~ m7r/(r2 — ri), where n and m are 
positive integers. We will consider only the lowest value of kr {m = 1) but allow n > 1. The 
initial perturbation is set to an approximate eigenmode appropriate for conducting boundary 
conditions: 

xu . . , n+r2-2r (r2-r)(r-ri) 

oBz — AsmkzZ oBj. — kzA cos kzZ- — oB^ — O 

r r 

5V, = ^cosA;^^ ^^^^^ 5Vr = KB sin k,z ~ ~ SV^ ^ 0. (10) 

r r 

Evidently, the fast-growing mode dominates the simulations no matter which n is used 
initially. Figure 4 compares the MRI growth rate obtained from the simulations with those 



-8- 



predicted by global linear analysis (Goodman & Ji 2002) as a function of magnetic Reynolds 
number. 

The radially global, vertically periodic linear analysis of Goodman & Ji (2002) found that 
the linear eigenmodes have boundary layers that are sensitive to the dissipation coefficients, 
but that the growth rates agree reasonably well with WKB estimates except near marginal 
stability. A comparison of the growth rates found by this analysis with those obtained from 
our simulations is given in Table 2. In the context of the simulations, "i?e = oo" means 
that the explicit viscosity parameter of the code was set to zero, but this does not guarantee 
inviscid behavior since there is generally some diffTision of angular momentum caused by 
finite grid resolution. Nevertheless, since the magnetic Reynolds number of the experiment 
will be about 20 and since Re/Rm ~ 10^, these entries of the table probably most closely 
approximate the degree of dissipation in the gallium experiment. In Table 2, the largest 
growth rate predicted by the hnear analysis has been marked with an asterisk (*). The 
simulations naturally tend to be dominated by the fastest numerical mode — that is, the 
fastest eigenmode of the finite-difference equations, which need not map smoothly into the 
continuum limit. Fortunately, as asserted by the Table, the fastest growth occurs at the 
same vertical harmonic n in the simulations as in the linear analysis. 

4. Nonlinecir Saturation 

As noted in §1, instabilities cannot easily modify the differential rotation of accretion 
disks because internal and magnetic energies are small compared to gravitational ones, and 
MRI is believed to saturate by turbulent reconnection (Fleming et al. 2000; Sano & Inutsuka 
2001). In Couettc fiow, however, the energetics do not preclude large changes in the rotation 
profile. As shown by Fig. 5), the differential rotation of the final state is reduced somewhat 
compared to the initial state in the interior of the flow, and steepened near the inner cylinder. 

4.1. Structure of the final state 

For moderate dissipation {Re,Rm < 10^), the final state is steady. Typical fiow and 
field patterns are shown in Figure 6. The poloidal flux and stream functions are deflned so 
that 

Vp = Vrer + V,e, = r-^e^X'V^, Bp = B,.e, + B,e, = r'^e^xV^ , (11) 

which imply V • Vp = and ^ • Bp = 0. [Our velocity field is slightly compressible, so 
that eq. (11) does not quite capture the full velocity field. Nevertheless, the error is small. 
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and $ is well defined by V^($e<^/r) = V X Vp with periodic boundary conditions in z and 
d'^jdz — ^ (ya. the cylinders.] 

The most striking feature is the outflowing "jet" centered near z = in Figure 6. The 
contrast in flow speed between the jet and its surroundings is shown more clearly in Figure 7. 
Figure 6 also shows that the horizontal magnetic field changes rapidly across the jet, which 
therefore approximates a current sheet. 

The radial flow speed in the jet scales with Rm as (Fig. 8), 

l^et « i?m-°-^^ (12) 
We find that the radial speed outside the jet scales similarly, 

External OC i^m-^'^^ ^ ^0.56_ (^g) 

Mass conservation demands that VjetWjet = Kxternai(2/i — VFjet), where Wjet is the effective 
width of the jet. Thus we can conclude that this width is independent of magnetic Reynolds 
number: 

W^jet OC i?m° (14) 

Additional support for this conclusion comes from the nearly equal scaling of Vj. and $ with 
Rm (Fig. 8), which indicates that the spatial scales in the velocity field are asymptotically 
independent of Rm. The toroidal fiow perturbation and toroidal field are comparable to the 
rotation speed and initial background field, respectively: 

1.18 < max ^ < 1.52, 0.28 < max < 0.56 (15) 

We emphasize that the scahngs (12)-(15) have been established for a limited range of flow 
parameters, 10^ < i?e, Rm < 10^-^. The jet is less well deflned at lower i?m, especially in the 
magnetic field. Extrapolation of these scalings to laboratory Reynolds numbers [Re > 10^) 
is risky, and indeed our simulations suggest that the final states are unsteady at high Re 
and/or high Rm (Fig. 11). 



4.2. Angular Momentum Transport 

Figure 9 displays the radial profiles of the advective, viscous, and magnetic torques 
integrated over cylinders coaxial with the boundaries: 



advective 



ive(r) = / 
J-h 



dz pr^VrV^ (16) 



-10- 



rmagnetic(?^) = \ dz \ ] (17) 

-h 
-h 

rvisnnnfi(^) / dz 



(18) 



3 9 (^V 

-r pf— 
i-h L or \ r 

rtotal('") — radvective('") ~l~ rmagnetic(^) ~l~ rviscous(^) (■^^) 

The advective and magnetic torques vanish at ri and r2 because of the boundary con- 
ditions but are important at intermediate radii. All components of the torque are positive 
except near r2- The total torque is constant with radius, as required in steady state, but 
increases from the initial to the final state (Figure 9). From Figure 10, we infer the scalings 

rfinal - Tinitial ^ ^^0.5^^0^ (20) 
i initial 

at least at Re, Rm > 10^. In fact, a better fit to the exponent of Re for Rm — 20 and 
Re > 10^ would be 0.68 rather than 0.5, but the exponent seems to decrease at the largest 
Re, and it is ^ 0.5 for Rm — 400, so we take the latter to be the correct asymptotic value. 

Representative runs are listed in Table 3. Additional runs have been carried out on 
coarser grids (smaller Nj., N^) to check that the values quoted for the torques are independent 
of spatial resolution to at least two significant figures in the laminar cases {Re, Rm < 10^) 
and to better than 10% in the unsteady cases where precise averages are difficult to obtain. 
In the latter cases, the quoted values in the last two columns have been averaged over radius 
but not over time. 



4.3. Interpretation of the final state 

The division of the fiow into a narrow outfiowing jet and a slower refiux resembles that 
found by Kagcyama ct al. (2004) in their hydrodynamic simulations [Fig. 3]. In that case, 
the jet bordered two Ekman cells driven by the top and bottom cndcaps. In the present 
case, however, Ekman circulation is not expected since the vertical boundaries are periodic, 
and we must look elsewhere for an explanation of the final state. 

Knobloch & Julien (2005, hereafter KJ) have proposed that axisymmetric MRI may 
saturate in a laminar flow whose properties depend upon the dissipation coefficients v h 7], 
with a large change in the mean rotation proflle, Vt{r). Although this mechanism of saturation 
probably cannot apply to thin disks for the reasons given in §1, it is consistent with some 
aspects of the final state of our Couette-fiow simulations: in particular, the scalings (12)-(13) 
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of the poloidal velocities with Rm; and the mean rotation profile does indeed undergo a large 
reduction in its mean shear, except near the boundaries (Fig. 5). 

One prominent difference between the final states envisaged by KJ and those found 
here is the axial lengthscale. KJ assumed the final state to have the same periodicity as the 
fastest-growing linear MRI mode, although they acknowledged that their theory does not 
require this. In our case, the linear and nonlinear lengthscales differ: whereas the fastest 
linear mode has three wavelengths over the length of the simulation (Table 2), the nonlinear 
state adopts the longest available periodicity length, namely that which is imposed by the 
vertical boundary conditions. Within that length, the flow is divided between the narrow 
jet and broad reflux regions. As discussed below, a third and even narrower reconnection 
region, whose width scales differently in Rm from that of the jet itself, exists within the 
jet. Another possibly important difference concerns the role of radial boundaries. KJ simply 
ignored these, yet our jet clearly originates at the inner cyhnder (Fig. 6). KJ's theory predicts 
that the poloidal flow should be proportional to i?e~^/^ as well as i?m~^/^ Yet, we flnd that 
V,,jct actually increases with Re, roughly as i^e"*"^/^, up to Re ~ 10^, above which it begins 
to decline and the flow becomes unsteady. 

The jet is probably the part of the flow that corresponds most closely to the "flngers" 
envisaged by KJ. Let us at least try to understand how the quantities in our jet scale with 
increasing Rm at fixed Re, even though it is more relevant to the experiment to increase Re 
at fixed and modest Rm (for the latter, see below). 

In steady state, the toroidal component of the electric field vanishes, — 0, because 
the flux through any circuit around the axis is constant. Consequently, 

The evidence from our simulations is that the peak values of $ and ^ scale as i]^^"^ and 77°, 
respectively, in the nonrestive limit rj ^ 0, Rm — > 00. The radial velocity Vr = r~^d^/dz 
also scales as r^^/^. In order that the two sides of eq. (21) balance, at least one of the 
derivatives of must become singular in the limit 77 — > 0. This appears to be the case. In 
fact, a comparison of the flux contours in Figures 6(a) and 12(a) suggests that a current 
sheet develops at the center of the jet. This is more obvious in the horizontal components 
of current density, Jr and J^, whose peak values we find to scale as oc ry-o-^^ pa Rm^^^ 
(Figure 13) and the maximum toroidal magnetic field near the current sheet scales as 

(X i?m°-^^ ^ Rm^'^ (22) 



From these scahngs one infers that the width of the current sheet scales as 77^/^. On the 
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other hand, the region defined by \Br\, \B^\ > \Bz\ appears to have a width oc 77°, hke that 
of the velocity jet. We call this the magnetic "finger" because of its form in Fig. 12. 

It is interesting to check whether these scahngs are consistent with the observation that 
the total torque (radial angular-momentum fiux) appears to be asymptotically independent 
of the resistivity. As — > 0, the advective torque oc J VrV^dz tends to zero since Vr oc 77^/^ 
and is presumably bounded by ~ rQi. The viscous contribution is always dominant 
near the cylinders but is reduced compared to the initial state at intermediate radii by the 
reduction in the vertically- averaged radial shear (Fig. 9). Since the total torque is larger 
in the final than in the initial state, a significant fraction of it must be magnetic, and this 
fraction should be approximately independent of 7] at sufficiently small r/. If ~ B^p oc rj^ 
within a vertical layer of width ~ 77^, the torque oc J B^B^dz oc rf^+v . Thus we expect 
y ~ —2x. In agreement with this, we have found that x ~ —1/2 and y ~ 1/3 in the current 
sheet, while in the finger, 2; ~ y ~ 0. 

One notices in Fig. 12(a)&(d) that the angular velocity is approximately constant along 
field lines — Q = ^(^) — as required by Ferraro's Law when the fiow is predominantly toroidal 
and the resitivity small. There must therefore be an outward centrifugal force along the lines 
in the magnetic finger, which in combination with the reconnection layer, presumably drives 
the residual radial outflow. Viscosity continues to be essential even as 77 — * because it 
is then the only mechanism for communicating angular momentum between field lines, and 
between the fluid and the cylinders; the distortion of the fleld enhances viscous transport by 
bringing into closer proximity lines with different angular velocity. 

To summarize, in the highly conducting limit Rm 00, Re =constant, there appear 
to be at least three main regions of the flow: (1) an "external" or "reflux" region in which 
the magnetic fleld is predominantly axial and the velocity predominantly toroidal, but with 
a small (oc 77^/^) radial inflow; (11) a "jet" or "flnger" of smaller but constant vertical width 
in which the flelds are mainly horizontal and there is a more rapid but still 0(77^/^) flow 
along fleld lines; (111) a resistive layer or current sheet at the center of the jet whose width 
decreases as 77^/^, across which the horizontal fields change sign. 

4.4. Simulations at small magnetic Prandtl number 

In the ongoing Princeton MRl experiment, the experiment material, liquid gallium, has 
kinematic viscosity v 2> x 10~^ cm^s~^ and resistivity 77 ~ 2 x 10^ cm^s~^. The typical 
dimensionless parameters are Rm ~ 20 and Re ~ 10^ at the dimensions and rotation speeds 
cited above. The magnetic Prandtl number Pr = Rm/Re fa 10"^ is very small. Rehable 
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simulations with Reynolds number as high as 10 are beyond any present-day computer, and 
small Pr presents additional challenges for some codes. 

Although our boundary conditions are not those of the experiment, we have carried 
out simulations at Rm — 20 and much higher Re in order to explore the changes in the 
flow due to these parameters alone. A simulation for Re — 25600 is shown in Figures 14 
& 15. All though this is still considerably more viscous than the experimental flow, it is 
clearly unsteady, like all of our simulations at Re > 3000. A narrow jet can still be observed 
in the poloidal velocities, but the poloidal field is only weakly perturbed at this low Rm: 

Since the Reynolds number of the experiment is much larger than that of our simulations, 
we can estimate the experimental torques only by extrapolation. Extrapolating according 
to eq. (20) from the highest--Re simulation in Table 3, one would estimate Ar/Finitiai ~ 35 
at Re - 10^. There are, however, reasons for caution in accepting this estimate. On the one 
hand, the experimental flow may be three-dimensional and turbulent, which might result in 
an even higher torque in the flnal state. On the other hand, the viscous torque in the initial 
state is likely to be higher than in these simulations because of residual Ekman circulation 
driven by the split endcaps. Nevertheless, we expect an easily measurable torque increase in 
the MRI-unstable regime. 

5. Conclusions 

In this paper, we have simulated the linear and nonlinear development of magnetoro- 
tational instability in a nonideal magnetohydrodynamic Taylor-Couette flow. The geometry 
mimics an experiment in preparation except in the vertical boundary conditions, which in 
these simulations are periodic in the vertical (axial) direction and perfectly conducting at 
the cylinders; these simplifications allow direct contact with previous linear studies. We have 
also restricted our study to smaller fiuid Reynolds number (Re), and extended it to larger 
magnetic Reynolds number (Rm), than in the experiment. We find that the time-explicit 
compressible MHD code ZEUS-2D, which is widely used by astrophysicists for supersonic 
ideal flows with free boundaries, can be adapted and applied successfully to Couette systems. 
MRI grows from small amplitudes at rates in good agreement with linear analyses under the 
same boundary conditions. Concerning the nonlinear flnal state that results from saturation 
of MRI, we draw the following conclusions: 

• Differential rotation is reduced except near boundaries, as predicted by Knobloch & 
Julien (2005). 
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• A steady poloidal circulation consisting of a narrow outflow (jet) and broad inflow 
is established. The width of the jet is almost independent of resistivity, but it does 
decrease with increasing Re. The radial speed of the jet oc Rm~^/'^. 

• There is a reconnection layer within the jet whose width appears to decrease oc RmT^^^. 

• The vertically integrated radial angular momentum flux depends upon viscosity but 
hardly upon resistivity, at least at higher Rm [eq. (20)]. 

• The final state is steady and laminar at Re, Rm < 10^ but unsteady at larger values 
of either parameter (Figs. 11 & 15.) 

• the final state contains horizontal fields comparable to the initial axial field for Rm > 
400, and about a tenth as large for experimentally more reahstic values, Rm ^ 20. 

We emphasize that these conclusions are based on axisymmetric simulations restricted 
to 10^ < Re, Rm < 10^'^, and that the boundary conditions are not realistic. This paper 
is intended as a preliminary exploration of MRI in the idealized Taylor- Couette geometry 
that has dominated previous linear analyses. We have not attempted to model many of the 
complexities of a realistic flow. In future papers, we will study vertical boundary condi- 
tions closer to those of the planned experiment; work in progress indicates that these may 
significantly modify the flow. 

The authors would like to thank James Stone for the advice on the ZEUS code. This 
work was supported by the US Department of Energy, NASA under grant ATP03-0084-0106 
and APRA04-0000-0152 and also by the National Science Foundation under grant AST- 
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Fig. 1. — Geometry of Taylor-Couette flow. In the Princeton MRI experiment, ri = 7.1 cm, 
r2 = 20.3 cm, h = 27.9 cm. 
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Fig. 2. — Radial profile of the azimuthal velocity for Re = 1. 
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Fig. 3. — Comparison with incompressible code at Re — 1600 : (a) Contours of toroidal 
velocity from Kageyama et al. (2004) (b) Results from ZEUS-2D with M = 1/4 




Fig. 4. — MRI growth rate versus Rm for conducting radial boundaries. Points: simulations. 
Curve: global linear analysis (Goodman & Ji 2002) with Re — 25, 600. 
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Fig. 5. — Angular velocity profile before and after saturation at several heights, for Re — 
Rm — 400. "Jet" is centered at 2; = (squares). 
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Fig. 6. — Contour plots of final-state velocities and fields. Re = 400, Rm = 400. (a) 
Poloidal flux function \1/ (Gauss cm^) (b) Poloidal stream function $ (cm^s~^) (c) toroidal 
field (Gauss) (d) angular velocity fl = r~^V^ (rads~^) 
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Fig. 7. — Radial velocity versus z for Re = 400, at several radii (cm): +, 8.42; *, 10.27; 
X, 11.98; A, 13.70; O, 16.87; □, 18.98. For clarity, only half the full vertical period (56 cm) 
is shown. Panel (a), Rm = 400; panel (b) Rm = 6400. 




Fig. 8. — Maximum radial speed in the jet (left panel) and maximum of poloidal stream 
function (right panel) vs. magnetic Reynolds number, for Re = 400. Powerlaw fits are shown 
as dashed lines with slopes —0.53 [left panel, eq. (12)] and —0.57 [right panel]. 
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Fig. 9. — 2;-integrated torques versus r. Re — 400, Rm — 400. Left panel: initial state; 
right: final state 
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Fig. 10. — Increase of total torque versus (a) Rm and (b) Re. In panel (b), dashed hues 
have slopes of 0.5 {Rm = 400) and 0.675 {Rm = 20). 
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Fig. 11. — Total toroidal magnetic energy vs. time at Re — 400. 
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Fig. 12. — Like Fig. 6, but for Rm = 6400, Re = 400. Symmetry about z = has not been 
enforced; the jet forms spontaneously ai z ^ ~20, but the whole pattern has been shifted 
vertically to ease comparison with Fig. 6. 
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Fig. 13. — Maximum radial current in the current sheet (left panel) and maximum of toroidal 
magnetic field (right panel) vs. magnetic Reynolds number, for Re — 400. Powerlaw fits are 
shown as dashed lines with slopes 0.46 [left panel] and 0.18 [right panel, eq (22)]. 
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Fig. 14.— Like Fig. 6, but for Re = 25600, Rm = 20. The flow is unsteady but closely 
resembles steady flows at lower Re for this Rm. 
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Fig. 15. — The z-averaged torques as in Fig. 9, but for the state shown in Fig. 14 {Re — 25600, 
Rm — 20). The radial variation of the total torque, though shght, testifies to the unsteadiness 
of the flow. 
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Table 1: Magnetic Diffusion Test 
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Table 2: Growth rates from semianalytic linear analysis vs. simulation. 
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Table 3: Increase of total torque versus Re and Rm. 
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